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Abstract We study the steady state of the abelian sandpile models with 
stochastic toppling rules. The particle addition operators commute with each 
other, but in general these operators need not be diagonalizablc. We use 
their abelian algebra to determine their eigenvalues, and the Jordan block 
structure. These are then used to determine the probability of different con- 
figurations in the steady state. We illustrate this procedure by explicitly 
determining the numerically exact steady state for a one dimensional exam- 
ple, for systems of size < 12, and also study the density profile in the steady 
state. 

Keywords Self-organized criticality, Stochastic Sandpile Model 



1 Introduction 

Sandpile models with stochastic toppling rules are important subclass of 
sandpile models [T]. The first such model was studied by Manna [2], and 
these are usually known as Manna models in the literature. They are able to 
describe the avalanche behavior seen experimentally in the piles of granular 
media much better than the deterministic models [Sj. Also, in numerical stud- 
ies, one gets better scaling collapse, and consequently, more reliable estimates 
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for the values of the critical exponents, than for models with deterministic 
toppling rules [HOa]. 

Unfortunately, at present, the theoretical understanding of models with 
stochastic toppling rules is much less than that of their deterministic counter- 
parts, e.g. the Bak-Tang -Wiesenfeld (BTW) model [6]. For example, there 
is no analogue of the burning test to distinguish the transient and the re- 
current states of a general Manna model. For the deterministic case, it is 
known that all the recurrent configurations occur with equal probability in 
the steady state. A similar characterization of the steady state is not known 
in the Manna case. The steady state has been explicitly determined only 
for the fully directed stochastic models [7UH] ■ In some cases, one can formally 
characterize the recurrent states of the model, e.g. the 1-dimcnsional Oslo rice 
pile model, but a straightforward direct depth-first calculation of the exact 
probabilities of different configurations in the steady state takes 0(exp(L 3 )) 
steps where L is the system length [9] . While the exact values of the criti- 
cal exponents have been conjectured for (1 + 1) dimensional directed Manna 
model |10pil|. the prototypical undirected Manna model in one dimension 
has resisted an exact solution so far p^2| J T3| J T4] . In higher dimensions, most 
of the studies are only numerical, and deal with the estimation of the critical 
exponents of avalanche distribution, and the universality class of the model 
[11II51II6U17UISU191I201EI] . 

While the original Manna model did not have the abelian property of 
the BTW model, one can construct stochastic toppling rules with abelian 
property [22]. In this paper, we discuss this abelian version of the stochastic 
Manna model. It is a special case of the more general Abelian Distributed 
Processors Model (ADP) pQ. We shall use the terms Deterministic Abelian 
Sandpile Models (DASM) and Stochastic Abelian Sandpile Models (SASM), 
if we need to distinguish between these two classes of models. 

We use the algebra of the addition operators to determine the steady state 
of the model. This algebraic approach provides a computationally efficient 
method to determine the Markov evolution matrix of the model. The addi- 
tion operators of SASM are not necessarily diagonalizable even if we restrict 
ourselves to the space of recurrent configurations. Using the abelian algebra 
we determine a generalized eigenvector basis in which the operators reduce 
to Jordan block form. We also define a transformation matrix between this 
basis and the configuration basis, and express the steady state in the latter 
basis. This procedure is illustrated by explicitly writing out the case of a one 
dimensional Manna model. In this special case, we can show that each Jordan 
block is at most of dimension 2. We determine the numerically exact steady 
state of the model for systems of size up to 12 and determine the asymptotic 
density profile by extrapolating the results. 

This paper is organized as follows: In section 2, we define the model 
precisely. In section 3, wc define the addition operators for the model, and 
discuss their algebra. Calculation of the eigenvalues and the Jordan block 
structure of the addition operators are given in section 4. The transformation 
matrix between the generalized eigenvector basis and the configuration basis 
is determined in section 5 and is used to determine the steady state vector 
in the configuration basis in section 6. The exact numerical determination of 
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the steady state is discussed in section 7 with some concluding remarks in 
section 8. 



2 The Model 

We define a generalized Manna model on a graph of N sites with a non- 
negative integer height variable Zi defined at each site i. Let the threshold 
height at i be zf, and the site is unstable if Zi > zf. If the system is stable, a 
sand grain is added at a randomly chosen site which increases the height by 
1. For each site i, there is a set of a™ ax lists E a ^ with a = 1, 2, • • • , u™ ax . 
If a site is unstable, it relaxes by the following toppling rule: we decrease its 
height by zf. Then, with probability p Q .i, we select the list E a ^, independent 
of any previous selections, and then add one grain to each site in that list. 
If a site occurs more than once in the list, we add that many grains to that 
site. 

Toppling at a site can make other sites unstable and they topple in their 
turn, until all the lattice sites are stable. It follows from the abelian property 
of the model that the probabilities of different final stable configurations are 
independent of the order in which different unstable sites are toppled. 

We illustrate these rules with some examples below. 

— Model A (The one dimensional Manna model): The graph is L sites on 
a line and zf = 2, for all sites. On toppling each grain is transfered to its 
neighbors with equal probability. Hence we have a™ ax = 3, for all i, with 

= {« i — 1) i — 1}, E 2 ,i = {i — 1, i + 1}, and E 3ji = {i + 1, i + 1} and 
Pi.i = P3.i = 1/4 an d P2.i = 1/2. Also grains can move out of the system 
if toppling occurs at a boundary site. 

— Model B (The one dimensional dissipative Manna model): Same as 
model A except that on toppling a grain can move out of the system with 
probability e. Then af lax = 6 and the lists of neighbors E\ = {i — 1, i — 1}, 
E 2 = {i - 1, i + 1}, E 3 = {i + 1, i + 1}, E 4 = {i- 1}, E 5 = {i + 1} and 
Eq = <P, where # is an empty set. The corresponding probabilities are 
Pi,i = P3,i = (1 - e) 2 / 4 , P2,i = (1 - e) 2 A P4,i = P5,i = e(l - e)/2 and 
P6,i = e 2 - 

In this case, one can use periodic boundary conditions, as there is dissi- 
pation at all sites, the steady state is critical only in the limit e — > 0. For 
the models A and B, it is easy to see that all stable configurations occur 
in the steady state with non-zero probability. We can also define stochas- 
tic models where the recurrent configurations form only an exponentially 
small fraction of all stable configurations. An example of this type is 

— Model C: The graph is a square lattice with N sites and zf = 2. Under 
toppling, with equal probability two particles are transfered to cither 
horizontal or vertical neighbors. Hence af lax = 2 with E\ t i = {i + e x , i — 
e^} and E 2 ,i = {i + e„, i - e y } with pi ti = p 2 ,i = 1/2. 

In the following we will mostly confine ourselves to Model A. The treat- 
ment of other cases presents no special difficulties. 
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3 The addition operators and their algebra 



N 



Let us denote the space of stable states as r spanned by Q = z\ basis 

vectors labeled by C. We define P(C,t) as the probability of finding the 
system in the basis C at time t. To each set {P(C, t)}, we associate a vector 
\P(t)) belonging to the vector space r, and write 



\P(t)) 



c 



P(C,t)\C). 



(1) 



Wc define the particle addition operators a^ for all i as linear operators acting 
on r as follows: Consider adding a sand grain at site i in a configuration C, 
and relaxing the system until a stable configuration is reached. For stochastic 
toppling rules, the resulting state is not necessarily a basis vector correspond- 
ing to a unique stable configuration, but a linear combination of them. If the 
resulting configuration is C with probability Pi(C'\C), we define 



ai \C) = Y,Pi{C'\C)\C), 



(2) 



for all C. Note that the action of any of these operators on a given configu- 
ration gives a unique probability state vector. 

Eq. ([2]) is a formal definition of the operators {a^}. One can think of these 
as J? x J? matrices, but it is quite non-trivial to actually determine the matrix 
elements Pi(G'\C) explicitly from the toppling rules. This is because of the 
non-zero probability of an arbitrary large number of toppling before a steady 
state is reached. 

For an example, consider the avalanches in model A for system of size L = 
3. Consider the 2 3 stable configurations as the basis vectors and denote them 
by their height values \z\,Z2,Zs). The action of a 2 on |0, 1,0) will generate 
a unstable state |0, 2, 0). Using the toppling rules we can write the following 
set of equations for three unstable states 



|0,2,0) 
|2,0,0) 
|0,0,2) 



j|2,0,0) 

1 

4 
1 

4 



||1,0,1) 



||0,0,2>, 



0,2,0) 
0,2,0) 



|0,1,0) 
|0, 1,0) 



0,0,0), 
0,0,0). 



(3) 



We see that there is a nonzero probability that the avalanche can continue 
for more than s toppling, for any finite s. e.g. in the sequence |0, 2, 0) — > 
|2,0,0) — > |0, 2, 0) • - - . Thus straight forward application of the relaxation 
rules do not result in a finite procedure to determine the unstable vector 
|0, 2, 0) in terms of the stable configurations. Instead, we have to write Eq. 
Q as a matrix equation 



(4) 





|0,2,0) 




|i,o,i) 


M 


|2,0,0) 




10,1,0) 




|0,0,2) 




10,0,0) 
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and then invert it. More generally, the determination of P(C'\C) involves 
working in a large space of unstable configurations. 

For example in model A, there are 2 L stable configurations, where each 
site has or 1 particle. Total number of particles is at most L. On adding 
one particle, the number of particles can become L + 1, where initially, only 
one site will have height 2. However, it is easy to verify that by toppling one 
can generate configurations where the number of particles at a site is much 
greater than 2. In fact, all the L + l particles could be at the same site. Then 
the total number of stable and unstable configurations Q 1 is the number of 
ways one can distribute L + l particles on L sites . It is easily seen that fl' 
varies as 4 L , and one needs to invert a matrix of size Oifl 1 x J?'). 

In this paper we will use the operator algebra to obtain an efficient method 
to determine the probabilities P(C'\C) explicitly which requires inverting a 
matrix only of size 2 L x 2 L . It has been shown [22] that the addition operators 
for different sites commute i.e. 



[a^ay] = 0, for all i,j. (5) 

Unlike the DASM, the inverse operators {a" 1 } for SASM need not exist, even 
if we restrict ourselves to the set of recurrent configurations. This is because 
among the recurrent states, one can have two different initial probability 
vectors that yield the same resultant vector. This makes the determination 
of the matrix form of the operators difficult for this model. 

Apart from the abelian property, the operators also satisfy a set of alge- 
braic equations. For simplicity of presentation, now on we consider zf = z c 
and p aj i = p a for all sites. Then consecutive addition of z c grains at a site 
ensures that the site will topple once and transfers z c grains to its neigh- 
bors, irrespective of the initial height. Then the operators obey the following 
equation 

af =^2p a a E «' i for 1 < i < N, (6) 

a 

where we have used the notation & E = Y\ a x for any list E, and 

xeE 

a, - 1, (7) 

for sites i outside the lattice. In particular for the examples in Section 2, 
these equations are as follows 

a, 2 = i(ai_i + a 4+1 ) 2 for Model A, (8) 

a 2 = [iy^a 4 _i + iy^a i+ i + el] 2 for Model B, and (9) 
a? = i(ai_ ex a i+ex + ai_e H ai+ e J for Model C. (10) 
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4 Jordan Block structure of the addition operators 

In general the matrices {a^} need not be diagonalizable. However, using the 
abclian property, we can construct a common set of generalized eigenvectors 
for all the operators {a,i} such that in this basis the matrices simultaneously 
reduce to Jordan block form. These generalized eigenvectors split the vector 
space r into disjoint subspaces, each corresponding to distinct set of eigen- 
values. There will be at least one common eigenvector in each subspace, for 
all the addition operators. 

Proof: Consider one of the operators, say ai . Let Ji be the subspace of 
r spanned by the (right) generalized eigen vectors of ai corresponding to 
the eigenvalue a\. There is at least one such generalized eigenvector, so Ji 
is non-null. We pick one of the other addition operators, say a2. From the 
fact that SL2 commutes with ai , it immediately follows that &2 acting on any 
vector in the subspace Ji leaves it within the same subspace. Diagonalizing 
&2 within this subspace, we construct a possibly smaller but still non-null 
subspace which is spanned by simultaneous eigenvectors of ai and a2 with 
eigenvalues a\ and 02. Repeating this argument with the other operators, one 
can construct vectors which are simultaneous eigenvectors of all the {a^}. 

Let |^>) be such an eigenvector, with 

a t \ip) = a^), for 1 < i < N. (11) 
Then from Eq.(6)the eigenvalues satisfy the following set of equations 

a\" =^p a a Ea ' i for 1 < i < N, (12) 

a 

where we have used the notation a E = Y\ a x, for any list E. 

xeE 

Rather than work with this general case, we will consider the special case 
in model A for simplicity. No extra complications occur in the more general 
case. Then, from Eq.Q, the corresponding eigenvalue equation is 

af = -(oi_i + a,+i) 2 , for 1 < i < L (13) 

These arc L coupled quadratic equations in L complex variables {ai}. We 
can reduce them to L linear equations by taking square root 

ViCti = -(a,-i + a i+i), ( 14 ) 

where r\i = ±1. The Eq. ([7]) sets the values for the eigenvalues of a and a £+1 
which are 

a 1 = a L+1 = l. (15) 

There are 2 L different choices for the set of L different 77's and for each such 
choice, we get a set of eigenvalues {ai}. In general, there will be degenerate 
sets of eigenvalues and the degeneracy arises if one of the ai is zero. Using 
the triangular inequality it is easy to show that 



2\ai\ < \cn-x\ + \a i+1 \, 



(16) 
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i.e. |aj| are convex functions of discrete variables i. Then, given the boundary 
condition in Eq. ()15|) . there could at most be one = in the solution for 
a given {?y;}, which means that each eigenvalue set can be at most doubly 
degenerate. 

Finding the number of degeneracies of solutions is interesting but difficult 
in general. We show that for L = 3 (mod 4) the number of such degenerate 
sets of eigenvalues > 2^- L+1 ^ 2 . 

Proof: Consider the system of length L = 4m + 3, with m being a non 
negative integer. For any given set {77;}, i = 1 to 2m + 2, it is possible to 
construct a solution {hi] of Eq. (fl4|) with i < 2m+2 which satisfies bo = 1 and 
&2m+2 = 0. Clearly, from Eq. (fl4j) . if we have the solution {a,i} corresponding 
to a particular set {77^ } , one can construct the solution {a^} corresponding to 
{r)j = —rjj} using = (— l) J 'oj. Using this symmetry we extend {bi} (i — 1 
to [L + l)/2) to form a set {a^} for i = 1 to L as follows: 

Oi = h for i < 2m + 2, (17) 

= (-l) 4 6 L+ i_, for i > 2m + 2. (18) 

This is a solution of Eq. (fT4")) for the set {77^} with 

/?■ = rii for i < 2m + 2, (19) 

= -r/i+i-i for i > 2m + 2. (20) 

and this solution {a^} satisfies the boundary conditions a = 1, ol+i = 1, and 
a2m+2 = (Fig.l). There are 2 2m+2 such solutions possible corresponding 
to all possible sets of {rj^}, and this gives the lower bound for the number of 
degenerate solutions. 

A direct numerical calculation for L < 20 shows that if L ^ 3 (mod 
4), all 2 L sets of eigenvalues are distinct. We present the degeneracies of 
the solutions in Table. 1. Calculation for simple choices of 77 shows that the 
degeneracies are possible only if L = 3 (mod 4). 
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32 


48 


288 


560 






19 


3024 
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288 





2432 



Table 1 Degeneracies arise if one of the a,i is zero in a solution of Eq. (|14|l , In 
the table, g denotes the total number of solutions with one of the at = i.e. the 
total number of degenerate sets of solution. Ni is the number of solutions with the 
eigenvalue a« = 0. Values for the other half of the system can be obtained using 
symmetry. 



For example, consider r)i = — 1 for i = L and for the rest of the sites it 
is 1. Then «j is of the form «j = 1 — od, for all i. If we want this to be zero 
for i = k, we must have a = 1/k. Then, requiring that the equation (fl~4|) be 
satisfied at i = L, gives 3L = 4k + 1, i.e. L = 3 (mod 4). Similarly for the 
set with r)L-i = —1 and 1 for rest of the sites imposes a condition on length 
3L = 8fc+ 1, which is also a subset of L = 3 (mod 4). Finding a general proof 
that degeneracies occur only if L = 3 (mod 4) remain an open problem. 

For each degenerate subspace there is a generalized eigenvector linearly in- 
dependent of the eigenvector corresponding to the eigenvalue of the subspace. 
In general, let us denote them by |{aj};n), where n = 1 for the eigenvector 
and n = 2 for the generalized eigenvector. For non-degenerate subspace n 
can only be 1. The vectors satisfy the following equations 

aj|{oj}; 1) = ai\{a,j}; 1), 

a i |{o J };2>=a i |{o J -};2) + a i |{o J };l) 1 (21) 

where a's are complex numbers. Then using the Eq. (|14|) it can be shown 
easily that a's satisfy the following equation 

This is similar to the Eq. (|14[) . except the boundary conditions which are 

a Q = a L+1 = 0. (23) 

For a given set of {?7i}, these are L simultaneous set of homogeneous linear 
equations which has infinitely many possible solutions. In order to get a 
single solution we choose ct!j = 1 if <Xi = 0, without loss of generality. This 
corresponds to choosing a particular normalization of the rank 2 eigenvectors. 
The solution of both the equations l|14p and <\'2'2\i can be easily obtained 
numerically. The generalized eigenvectors and the Jordan block form of the 
addition operators for the system of size L = 3 are given in the appendix. 

5 Matrix representation in the configuration basis 

Given the well defined action of the addition operators on the generalized 
eigenvectors it is possible to define a transformation matrix M between the 
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configuration basis and the generalized eigenvector basis. 

|{*}>=E M {*}.J#i>> (24) 
o 

where |{^}) is the basis vector of r corresponding to the height configu- 
ration {zi} and \ipj) is the jth generalized eigenvector. Let us express the 
configuration |{0}}, with all sites empty, as a linear combination of all the 
generalized eigenvectors. 

|{0}> (25) 

3 

where CjS are constants. Then all the stable configurations can be obtained 
by adding grains at properly chosen sites in |{0}). 

\{zi}) = n < ! k°» = e °> n < ! w> w 

i j i 

and hence 

M {a4}J =({« i }in a ?i^>- ( 27 ) 

i 

The action of the addition operators on the generalized eigenvectors, for 
example Eq. pTj) for model A, would generate the elements of the matrix M. 
Given M, we can get the eigenvectors of a^, in the configuration basis, in 
particular, the steady state vector, by the inverse transformation 

\i> j )=M- 1 \{z i }}. (28) 

The addition operators in the configuration basis are obtained using the 
similarity transformation Ma/M _1 . An explicit form of M for model A of 
length L = 3 is given in the appendix. 



6 Determination of the steady state vector 

The time-evolution of the system is Markovian and the evolution operator 
W is defined by the master equation 

\P(t+l))=W\P(t)), (29) 

where \P(t)) and \P(t + 1)) are the state of the system at time t and t + 
1, respectively. We can write the time-evolution operator in terms of the 
addition operators as 

W=i^>,. (30) 

i 

Then the common eigenvector of all the addition operators corresponding to 
eigenvalue 1 is the steady state vector of the system. The steady state vector 
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Fig. 2 The ratio of the probability of the most probable configuration Cmax (all 
occupied) and the least probable configuration C m in (all sites empty) plotted as a 
function of the system length L. The fitting function f(x) = a — bL + cL log L, with 
a = 1.50, b = 0.80 and c = 0.94. 



can be determined in the stable configuration basis using the matrix M . 
For model A of length L = 3 the steady state vector is 

13 1 47 3 

|5>= ^|0,0,0) + ^|1,0,0) + -|0,1,0) + -|1,1,0) 

+ i|0,0 > l) + H|i j0> l)+l[0, 1,1) + ^11,1,1), (31) 

where the stable configurations are denoted by \zi,Z2,Z3) with Zj as the 
height of the ith site. The amplitude of each term in the expansion are the 
probability of finding the steady state in the corresponding height configu- 
ration. 



7 Numerical Results 

Here we numerically calculate the exact steady state of model A for different 
system length and discuss its properties. As shown in Eq. fTT)) and Eq. 
the eigenvalues {a{\ and the off-diagonal matrix elements {c^} form sets of 
linear equations for a given set of {r)i}. We solve them by LU decomposition 
method. Because of the tridiagonal structure of the equations, only 0(2 L ) 
number of steps are required to get the solution. The maximum number of 
steps (0(2 3L )) are required for the inversion of the transformation matrix 
M. We have used the Gauss-Jordan elimination method for the inversion. It 
is important to note that, the maximum system length L, possible to treat 
by this method, is determined by the limited memory size of the computers, 
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and not by the computation time. Using desktop computers we were able to 
determine M exactly for systems of size L < 12. 

We note that as L is increased, the second largest eigenvalue of W tends 
to 1/2. Thus, the gap between the largest and the next largest eigenvalue of 
the relaxation matrix does not tend to zero. This gap measures the relaxation 
time of the system in terms of the macro-time unit of interval between ad- 
dition of grains. However, the average duration of an avalanche measured in 
terms of micro-time unit of duration of a single toppling event does diverge, 
as system size increases. 

An interesting question is the extent of variation between probabilities 
of different configurations in the steady state. In the one-dimensional Oslo 
model, for a system of L sites, the ratio of probabilities of the most probable 
to the least probable configuration varies as exp(L 3 )[9]. However in model A, 
we find that the ratio is not quit as large, and it only varies approximately 
as exp(0.94ilogi) (FigEJ) for large L. 

This suggests that possibly the exact steady state has a product measure. 
To check this we define a product basis = TT |V>i), where could be 

i 

any one of the two orthogonal vectors 

|1') = cos&|l) + sin&|0>, 

|0') =sin&|l)-cos&|0), (32) 
with <f>i a real number. Then in this basis the steady state can be written as 

\S)=Y.P^')W)- (33) 

We choose {4>i} so that the ratio between the amplitudes of basis vectors 
with next-largest and largest amplitudes becomes as small as possible (this 
would become zero, if the state was a product measure state). In Figj3l 
we have plotted for system of size L = 12, the relative amplitudes in both 
configuration basis and the optimized product basis as a function of the rank 
of the basis vectors with the vectors arranged in decreasing orders of their 
amplitudes. In the optimized basis the second highest probability is only 10 
times smaller than the highest probability. This implies that the steady state 
measure is not a product measure. 

The steady state density for different sites are plotted in Fig.|4]for different 
system sizes. Amongst the different fitting forms that we tried, the following 
functional form gives the best fit 

+ , , \ : ^ (34) 



Pl(x) poo (x + d) u ^- (L + l- x + d)'^- 

where poo, b, v± and d are real numbers. Using this functional form the steady 
state particle density averaged over all sites for system of size L can be written 
as 

Tl = T^ + JlTW 1, (35) 
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Fig. 3 The amplitudes, normalized with its largest value, corresponding to the 
basis vectors in the steady state plotted as a function of the rank of the basis 
vectors. The vectors are arranged in decreasing orders of their amplitudes. The 
plot is given for the configuration basis and the optimized basis for model A of size 
L = 12. 
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Fig. 4 Average steady state density at site i for model A of different length 

L. 



where B is a real number and poo is the asymptotic value of the average 
particle density. The exact value of pL and the particle density at the central 
site Pi(cc m ) are listed in the Table 2 for different system sizes. The sequential 
fitting method is used to find the values of poo, B, v± and S from these data. 
For a given choice of S, these values are obtained numerically by solving the 
Eq. (|35|) for three consecutive lengths L — 1, L and L + 1. Best convergence 
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L 


PL 




2 


0.583333 


0.583333 


3 


0.634354 


0.709184 


4 


0.669262 


0.737000 


5 


0.695210 


0.769704 


6 


0.715472 


0.786491 


7 


0.731879 


0.805897 


8 


0.745514 


0.816009 


9 


0.757080 


0.827217 


10 


0.767051 


0.834600 


11 


0.775760 


0.842665 


12 


0.783451 


0.848054 



Table 2 The values of particle density in the steady state for the model A of 
different length L. Here pL denotes the steady state particle density averaged over 
all sites and pL(x m ) denotes the steady state particle density at the central site. 



L 


1/Poo 


B 


v±. 


3 


1.061 


1.128 


0.656 


4 


1.049 


1.132 


0.641 


5 


1.053 


1.132 


0.646 


6 


1.049 


1.131 


0.640 


7 


1.050 


1.131 


0.641 


8 


1.049 


1.130 


0.639 


9 


1.049 


1.130 


0.639 


10 


1.049 


1.130 


0.639 


11 


1.049 


1.130 


0.639 



Table 3 The sequential fit of the functional form in Eq. (|35|l to the data for average 
particle density for model A of different length L given in Table 2. 



of the values of p o 1 B and v± are obtained for 5 = 1.1, which are tabulated 
in Table 3. The asymptotic value of the average particle density converges 
to poo = 0.953 which is close to the more precise estimate 0.94885(7), from 
Monte Carlo simulations [12] . 



8 Concluding remarks. 

For a general SASM with N sites, the calculation of eigenvalues involves 
solving N coupled polynomial equations in N variables. This can be done 
in polynomial time in J?, the number of stable configurations of the model. 
These are then used to construct the transformation matrix M of size f2y.il. 
Finally inverting the matrix M gives us the eigenvectors of the evolution 
operator, in particular the steady state. 

Of course, to determine the steady state of any Markov chain on Q states, 
we need to determine the eigenvectors of the evolution matrix of size fl x fl. 
The point here is that the specification of the toppling rules does not directly 
specify the evolution matrix, and determining the matrix elements of the 
latter from the toppling rules is computationally very nontrivial. Using the 
abelian property, we are able to tackle this problem. 
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For a generic model with some parameters, e.g. the model B, except for 
special symmetries, one does not expect degeneracies in eigenvalues to occur 
for a generic value of the parameters. For special values of the parameters, 
if there is a non-trivial Jordan block structure of the evolution operator, it 
would show up in the timc-dcpcndcnt correlation functions of the model by 
the presence of terms of the type iexp(— Xjt), in addition to the usual sum 
of terms of the type exp(— Xjt). 

In particular we have explicitly calculated the steady state for a specific 
model (model A in section 2) of length L< 12. Extrapolating the results we 
determined the asymptotic density profile in the steady state. The power-law 
profile of deviations from the mean value near the ends would be important 
for determining the avalanche exponents of the model |23j . This remains an 
interesting open problem. 

Acknowledgements The work of DD is supported in part by a J. C. Bose Fel- 
lowship of the Government of India. 



A Appendix 

Here we give some details of the explicit calculation of the steady state, and the 
matrix representation of addition operators for model A of length L = 3. 

The eight sets of eigenvalues obtained by solving Eq, (jl4[) are (1, 1, 1), ( — 1, 1, — 1), 
(|, — |, |), (— |, — |, — |), (|, 0, — |), and ( — |, 0, |) with the last two sets repeated 
twice. 

For writing the matrix structure of the addition operators, we choose the order 
of the eigenvectors same as the order of the eigenvalues mentioned ab ove. For the 
degenerate subspace we order the eigenvector |{hi}; 1), defined in Eq. (|21|l . before 
the generalized eigenvector |{ai}; 2). Then in this basis the matrices corresponding 
to the addition operators have the following Jordan block form. 
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(38) 



The transformation matrix M, discussed in section 5, between the generalized eigen- 
vector basis and the configuration basis has the following form 



M : 
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1 
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1/3 
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1/2 


-1 




i 


1 


-1/3 


-1/3 
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-1 


-1/9 


1/9 


1/2 
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-1/3 
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-1/27 


-1/27 


-1/4 





-1/4 


o ) 



(39) 



where the configuration basis vectors are chosen in the following order (0,0,0), 
(1,0,0), (0,1,0), (1,1,0), (0,0,1), (1,0,1), (0,1,1), and (1,1,1). The matrix is 
non-singular, and the inverse can be calculated numerically. Using the similarity 
transformation Ma/M -1 we find matrix representation of the addition operator 
ai in the configuration basis. 
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The other operators can also be determined similarly. 
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